Glacier mass balance
We simplify the glacier mass balance to the following equation
\[ \Delta S = P-M = A_{imbal} \]
where the change of water storage (\(\Delta S\)) is equal to the precipitation (\(P\), assumed solid) minus the glacier melt (\(M\)). Typically melt excesses precipitation and we have negative imbalance ablation, indicating glacier storage loss. The glacier mass balance is calculated in annual time steps. We thereby refer to the hydrological year starting on October 1st of the previous year to take a full accumulation and ablation season into account.
Temperature index model
Hock (2003) describes several variants of the temperature index model for simulating glacier melt. The package implements the temperature index melt model described in (Hock 2003) in the function .
\[ M = \biggl\{ \begin{array}{l, l} 0, & T < T_{threshold} \\ f_{M} \cdot \left( T - T_{threshold} \right), & T >= T_{threshold} \end{array} \]
where \(M\) is the glacier melt in \(mm/d\), \(T\) is the daily average temperature in \(^{\circ} C\). The two parameters \(f_{M}\) and \(T_{\text{threshold}}\) refer to the melt factor and the threshold temperature above which glacier melt occurs and need to be calibrated. They have the units \(\frac{mm}{^{\circ} C \cdot d}\) and \(^{\circ} C\) respectively. Glacier melt is calculated in daily time steps.
Glacier volume development
As glaciers melt, their volume changes. This has to be taken into account for the long-term simulation of glacier discharge. To determine the initial glacier volume we use the area of the geometry of the RGI v6.0 data set and multiply it with the average thickness of the glacier. Glaciers larger than 1km2 we sub-divide into elevation bands of 100m altitude to account for elevation dependent temperature forcing. The large glaciers melt from the lowest elevation band to the highest elevation band whereby the glacier melt is subtracted from the glacier volume of lowest elevation band which is still glacierized. The small glaciers are not spatially discretized and thus they are melted homogeneously.
For annual time step \(t\), the evolution of the glacier volume is calculated as follows:
\[ A(t) = \text{glacierArea_RGIF}\bigl(V(t)\bigr) \]
\[ Q_{glacier}(t) = M(t) \cdot A(t) \]
\[ V(t+1) = V(t) + \Delta S = V(t) + \text{glacierImbalAbl}\bigl(M(t)\bigr) \]
Note that the imbalance ablation (or \(\Delta S\)) is negative for glacier mass loss and positive for glacier accumulation.
\(Q_{\text{glacier}}\) can be calibrated against glacier discharge derived from the Miles & Hugonnet data sets. The automated calibration is not yet included in the package.
The function is an empirical scaling function analogue the inverse of the scaling function derived by Erasov (1968) but based on the modern RGI v6.0 glacier geometries and the glacier thickness data set by Farinotti et al. (2019). The package implements volume-area and area-volume scaling functions based on both, Erashov and RGI-Farinotti data, allowing the estimation of glacier areas based on glacier volumes ( and ) and estimations of glacier volumes based on glacier areas ( and ).
The function is an empirical scaling function relating glacier imbalance ablation, i.e. the glacier storage loss glacier melt. It is derived from the glacier discharge data set by Miles et al. (2021) and the glacier thinning data set by Hugonnet et al. (2021).
Propagation of uncertainty
The initial glacier volume of each glacier is attributed an uncertainty of p/m 26%. This number is based on the average uncertainty of the glacier volume per RGI region reported in Farinotti et al. (2019) (Table 1). The error is the same for RGI regions 13 and 14. The uncertainty of the area of the RGI v6.0 glacier outlines is not known. We therefore assume that the error of the glacier volume stems to 50% from the estimation of the glacier thickness and to 50% from the glacier area. We further assume that the errors of the glacier area and glacier thickness are un-correlated and can thus estimate the uncertainties of the glacier area data and the glacier thickness data to be p/m 13% each.
For the non-linear relationships and , the standard deviation of the residuals of the fit was computed. The estimated error of the fit is assumed to be equal plus/minus twice the standard deviation of the residuals and yields 31% and 53% respectively. The residuals are not normally distributed and their actual distribution in unknown. Nevertheless, we apply the principles of linear error propagation for the sake of simplicity. We assume uncorrelated errors and therefore add the error of the function input to the error of the fit to obtain the error of the function output.
\[ \varepsilon_{V} = \varepsilon_{A} + \varepsilon_{\text{glacierVolume_RGIF}} = 0.26 + 0.31 = 0.57 \]
\[ \varepsilon_{A} = \varepsilon_{V} + \varepsilon_{\text{glacierArea_RGIF}} = 0.26 + 0.53 = 0.79 \]
We do not have an estimate for the error of the temperature index model. We assume a conservative relative error of 2, indicating that the estimated glacier melt is within a range of plus/minus 2 times it’s value.
The error of the imbalance ablation amounts to 73%. Assuming independent errors from the temperature index model for glacier melt and the scaling function between imbalance ablation and glacier melt, the error of the imbalance ablation amounts to approximately:
\[ \varepsilon_{Q_{\text{imb,melt}}} = \varepsilon_{Q_{\text{M}}} + \varepsilon_{\text{glacierImbalAbl}} = 2 + 0.73 = 2.73 \]
The errors stated above relate to initial estimates of the glacier areas, volumes, total ablation and imbalance ablation. The uncertainties are not propagated through each time step of the model as the errors depend on the errors of the previous time step, i.e. they are not uncorrelated over time. In any case, the uncertainty margins for any glacier melt modelling done with the presented method are large.
# Relative error estimates for the initial state in %
error_stats = tibble::tibble(
sV = 0.26, # Glacier volume
sA = 0.13, # Glacier area
sTh = 0.13, # Glacier thickness
sMelt = 2, # Glacier melt
sVolumeRGIF = 0.31, # Calculating glacier volume from glacier area
sAreaRGIF = 0.53, # Calculating glacier area from glacier volume
sImbal = 0.73, # Calculating glacier imbalance ablation from glacier melt
sQglgrowth = 1.6, # Calculating total ablation from glacier thinning
sQglmelt = 1.6) Demonstration
We demonstrate the above described method with the data from the Atbashy basin, a tributary of the Naryn river. As for the previous chapter you can download all required data from this link.
library(tmap)
library(timetk)
library(sf)
library(raster)
library(exactextractr)
library(tidyverse)
library(lubridate)
devtools::install_github("hydrosolutions/riversCentralAsia")
library(riversCentralAsia)
# Stochastic optimisation algorithm for model calibration
library(GA)
# Path to the data directory downloaded from the download link provided above.
# Here the data is extracted to a folder called atbashy_glacier_demo_data
data_path <- "../../atbashy_glacier_demo_data/"Spatial discretization
Glacier melt is calculated on elevation bands. The following code chunks show how to derive elevation bands on glaciers.
TODO!
Below figure shows the division of the larger glaciers into elevation bands. Glaciers smaller than 1 km2 are not sub-divided to limit the use of computational resources. The area threshold of 1 km2 is arbitrary but similar area thresholds have been used in the literature (e.g. Miles et al. (2021)).
# Load data
## Digital elevation model (DEM)
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif"))
## Glacier outlines from the Randolph Glacier Inventory (RGI) v6.0
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
## Outlines of hydrological response units for the modelling of glacier discharge.
rgi_elbands <- st_read(paste0(data_path,
"GIS/rgi_glaciers_atbaschy_el_bands.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
# Load a pre-processed raster file with glacier thickness in the Atbashy basin
# (see vignette [glaciers 01]{glaciers-01-intro.html} for details on the pre
# -processing)
glacier_thickness <- raster(
paste0(data_path,
"GLACIERS/Farinotti/pre-processed_glacier_thickness.tif"))
# Only keep the variables we need for this analysis
rgi_elbands <- rgi_elbands |>
dplyr::select(RGIId, Area, elvtn_b) |>
rename(Area_tot_glacier_km2 = Area) |>
mutate(ID = paste0(RGIId, "_", elvtn_b))
# Get mean elevation of each glacier/elevation band from DEM
rgi_elbands$z_masl <- exact_extract(dem, rgi_elbands, "mean", progress = FALSE)
# Update the glacier area within the basin boundaries
rgi_elbands$A_km2 <- as.numeric(st_area(rgi_elbands))*10^(-6)
glaciers <- unique(rgi_elbands$RGIId)
# Calculate the average glacier thickness for each elevation band.
rgi_elbands$thickness_m = exact_extract(glacier_thickness,
rgi_elbands, "mean", progress = FALSE)
rgi_elbands$thickness_m[is.na(rgi_elbands$thickness_m)] <- 0
tmap_mode("view")
tm_shape(dem, name = "DEM") +
tm_raster(n = 6,
palette = terrain.colors(6),
alpha = 0.4,
legend.show = TRUE,
title = "Elevation DEM (masl)") +
tm_shape(rgi, name = "RGI v6.0") +
tm_borders(col = "gray", lwd = 0.4) +
tm_shape(rgi_elbands, name = "Elevation bands") +
tm_polygons(col = "elvtn_b", n = length(unique(rgi_elbands$elvtn_b)),
palette = "Blues",
lwd = 0.2, alpha = 0.8, title = "Elevation band ID") +
tm_scale_bar(position = c("right", "bottom")) Observed glacier melt
Observed average glacier thinning rates by Hugonnet et al. (2021) are shown in the figure below. 30 out of 135 glaciers in the Atbashy basin show increase of glacier thickness between 2000 and 2020.
# Load the glacier thinning data set and filter it to the glaciers in the
# catchment of interest.
hugonnet <- read_csv(paste0(
data_path, "/GLACIERS/Hugonnet/dh_13_rgi60_pergla_rates.csv"),
show_col_types = FALSE)
hugonnet <- hugonnet |>
dplyr::filter(rgiid %in% rgi$RGIId) |>
tidyr::separate(period, c("start", "end"), sep = "_") |>
mutate(start = as_date(start, format = "%Y-%m-%d"),
end = as_date(end, format = "%Y-%m-%d"),
period = round(as.numeric(end - start, units = "days")/366))
glaciers_hugonnet <- rgi |>
left_join(hugonnet |> dplyr::select(rgiid, area, start, end, dhdt, err_dhdt,
dvoldt, err_dvoldt, dmdt, err_dmdt,
dmdtda, err_dmdtda, period),
by = c("RGIId" = "rgiid"))
# Add glacier discharge by Miles et al.
miles <- read_delim(paste0(data_path,
"GLACIERS/Miles/Miles2021_Glaciers_summarytable_20210721.csv"),
show_col_types = FALSE, delim = ",") |>
dplyr::filter(RGIID %in% rgi$RGIId & VALID == 1)
glaciers_hugonnet <- glaciers_hugonnet |>
left_join(miles |> dplyr::select(RGIID, totAbl, totAblsig, imbalAbl,
imbalAblsig),
by = c("RGIId" = "RGIID")) |>
mutate(Qgl_m3a = ifelse(is.na(totAbl), glacierDischarge_HM(dhdt), -totAbl))
observed_melt <- glaciers_hugonnet |>
st_drop_geometry() |>
dplyr::filter(period == 1) |>
dplyr::select(RGIId, start, dhdt, Qgl_m3a, Area) |>
mutate(Area_m2 = Area * 10^6) |>
dplyr::select(-Area) |>
distinct()
totQgl <- glaciers_hugonnet |>
st_drop_geometry() |>
dplyr::filter(period == 20) |>
mutate(Qgl_m3a = ifelse(Qgl_m3a < 0, 0, Qgl_m3a),
dVweqda_m3a = ifelse(dmdtda<0, -dmdtda, 0) * area) |>
summarise(Qgl_m3a = sum(Qgl_m3a),
dVweqda_m3a = sum(dVweqda_m3a)) |>
transmute(Qgl_mcma = Qgl_m3a * 10^(-6),
Qgl_m3s = Qgl_m3a / (60*60*24*365),
dVweqda_mcma = dVweqda_m3a * 10^(-6),
dVweqda_m3s = dVweqda_m3a / (60*60*24*365))
tmap_mode("view")
tm_shape(dem, name = "DEM") +
tm_raster(n = 6,
palette = terrain.colors(6),
alpha = 0.8,
legend.show = TRUE,
title = "Elevation (masl)") +
tm_shape(glaciers_hugonnet |> dplyr::filter(period == 20) |>
mutate(Qgl_mcma = Qgl_m3a * 10^(-6)),
name = "Glacier discharge\n[10^6 m3/a]") +
tm_polygons(col = "Qgl_mcma", lwd = 0.4,
title = "Qgl (10^6 m3/a)") +
tm_scale_bar(position = c("right", "bottom"))The average annual contribution of glacier melt to river runoff is 1.9704266 m3/s. The annual norm discharge at Atbashy gauging station is 16.6 m3/s which yields an average annual glacier contribution of about 11.8700399 % which is reasonable given the high median altitude of the basin of 3470 masl.
Forcing
There will be a separate Vignette to demonstrate how to prepare the forcing. The methodology is described in the CAHAM book. For now, we provide you with the pre-processed forcing data and give you just a brief overview over the data source. We have sub-divided the glaciers into elevation bands to be able to reproduce the elevation-dependent temperature forcing on the glacier.
Historical temperatures
As meteorological data for high elevations in Central Asia is very scarce we use the CHELSA v2.1 data set (Karger et al. 2017). This is a global data set of forcing data for hydrolgical models based on ERA5 but corrected for biases and especially suited for high elevations. The daily CHELSA forcing has been cut to the Central Region by the originator of the data set, D. Karger, WSL and extracted to the hydrological response units of the glaciers in the Atbashy basin by T. Siegfried.
hist_obs <- readRDS(file = paste0(data_path, "CLIMATE/hist_obs_glacier_tas.rds"))
# Plot the temperature time series for a given glacier/elevation band
hist_obs[, c(1, 2, 5, 10)] |>
pivot_longer(-date, names_to = "layer", values_to = "Temp_degC") |>
mutate(layer = as.numeric(gsub(pattern = "_", "",
str_extract(layer, "_\\d$")))) |>
plot_time_series(
.date_var = date,
.value = Temp_degC,
.color_var = layer,
.color_lab = "El. band",
.plotly_slider = TRUE,
.smooth = FALSE,
.line_size = 0.5
)Daily temperature time series of the lowest, a middle-elevation and the highest elevation band of glacier RGI60-13.08930 in the Atbashy basin. Data source: CHELSA v2.1.
Future temperatures
Future temperature development per glacier or elevation band is extracted from the 3 CMIP6 GCM models with highest priorities for the region (TODO Add names and citations) downloaded from COPERNICUS. We take 4 socioeconomic scenarios into account, covering 4 different emission scenarios. The temperatures of the climate models are bias corrected using the CHELSA data and a quantile mapping method. The figure below shows the development of the temperature in the lowest elevation band of glacier RGI60-13.08930 as an example.
fut_sim <- readRDS(file = paste0(data_path, "CLIMATE/fut_sim_glacier_tas_qmapped.RDS"))
# Plot the temperature time series for a given glacier/elevation band
glacier <- "RGI60-13.08930_1"
# Extract the temperature for the selected glaciers for all GCMs and SSPs
scenarios <- names(fut_sim)
fut_temp <- NULL
for (scenario in scenarios) {
fut_temp <- rbind(fut_temp,
fut_sim[[scenario]] |>
dplyr::select(Date, all_of(glacier)) |>
mutate(Scenario = scenario))
}
fut_temp <- fut_temp |>
mutate(Hyear = hyear(Date)) |>
group_by(Hyear, Scenario) |>
summarise(Date = first(Date),
Temp = mean(get(glacier))) |>
separate(Scenario, into = c("GCM", "SSP"), sep = "_") |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear),
GCM != "IPSL-CM6A-LR")
# Plot annual data for better readability
ggplot() +
geom_line(data = hist_obs |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear) |>
summarise(date = first(date),
Temp = mean(get(glacier))) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)),
aes(date, Temp)) +
geom_line(data = fut_temp, aes(Date, Temp, colour = GCM,
linetype = SSP)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "T [deg C]", title = glacier) +
theme_bw()
Annual historical (CHELSA) and future (CMIP6, Copernicus) temperature for the lowest elevation band of glacier RGI60-13.08930.
Calculate daily glacier melt
The code snippet below illustrates how to calculate the daily glacier melt in a river catchment with default parameters for the temperature index model. The figure shows daily melt rates for the different elevation bands of one of the larger glaciers (RGI60-13.08930). The highest melt occurs at low elevation bands. Glacier melt only occurs during the summer months when the average daily temperature is above the threshold temperature for glacier melt.
# Default parameters of the temperature index model for glacier melt.
MF_small = 2 # Melt coefficient for small glaciers (smaller than 1km2)
MF_large = 0.1 # Melt coefficient for large glaciers (larger or equal than 1km2)
threshold_temperature = 0 # Default temperature threshold for glacier melt
Area <- rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, A_km2) |>
pivot_wider(names_from = ID, values_from = A_km2)
# Assign different melt factors to large and small glaciers.
MF <- Area |>
mutate(across(everything(), ~MF_large),
across(ends_with("_1"), ~MF_small))
melt <- glacierMelt_TI(temperature = hist_obs |> dplyr::select(-date),
MF = MF,
threshold_temperature = threshold_temperature)
melt <- as_tibble(melt) |>
mutate(date = hist_obs$date) |>
relocate(date, .before = where(is.numeric))
ggplot(melt[, 1:9] |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
left_join(rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, z_masl), by = "ID") |>
mutate("Mean\nelevation\n[masl]" = as.factor(round(z_masl)))) +
geom_line(aes(date, Melt, colour = `Mean\nelevation\n[masl]`)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "Melt [mm/d]") +
theme_bw()
Daily glacier melt per elevation band.
Compare to measured glacier melt
Aggregate the daily glacier melt to annual data and compare it to the observed glacier melt.
melt_a_eb <- melt |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt))|>
ungroup()
melt_a <- melt_a_eb |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(Hyear, RGIId) |>
summarise(date = as_date(first(date)),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt)) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |>
left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(),
by = "RGIId") |>
# Assume constant glacier area over the comparison period.
mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6,
.lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),
.lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a),
.ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))
melt_obs_a <- observed_melt |>
mutate(.lb_Qgl_m3a = ifelse(
dhdt > 0,
ifelse(Qgl_m3a*(1-error_stats$sQglgrowth)<0, 0,
Qgl_m3a*(1-error_stats$sQglgrowth)),
ifelse(Qgl_m3a*(1-error_stats$sQglmelt)<0, 0,
Qgl_m3a*(1-error_stats$sQglmelt))),
.ub_Qgl_m3a = ifelse(dhdt > 0,
Qgl_m3a*(1+error_stats$sQglgrowth),
Qgl_m3a*(1+error_stats$sQglmelt)))
melt_obs_a <- melt_obs_a |>
dplyr::filter(RGIId %in% glaciers[c(6, 7)])
ggplot() +
geom_ribbon(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6,7)]),
aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a,
fill = RGIId), colour = NA, alpha = 0.1) +
geom_ribbon(data = melt_obs_a,
aes(start, Qgl_m3a,
ymin = .lb_Qgl_m3a,
ymax = .ub_Qgl_m3a,
colour = RGIId, linetype = "obs", fill = RGIId),
size = 0.2, alpha = 0.1) +
geom_line(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6,7)]),
aes(date, Qgl_m3a, colour = RGIId, linetype = "sim")) +
geom_line(data = melt_obs_a, aes(start, Qgl_m3a, colour = RGIId,
linetype = "obs")) +
labs(x = "Date", y = "Melt [m weq/a]") +
scale_linetype_manual(name = "Source",
values = c("sim" = 1, "obs" = 2)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ggtitle(paste0("MF: ", MF, "Tth: ", threshold_temperature)) +
theme_bw()
Daily glacier melt per elevation band.
The temperature index model can be calibrated with the specific glacier volume change provided by Hugonnet et al. (2021) (see \(\ref{fig:glaciermeltcomparison}\)).
Calibration of the temperature index model
The 2 parameters of the temperature index model can be calibrated using observed glacier melt rates. As mentioned above, the observed glacier melt rates are simulated total ablation by Miles et al. (2021) where available or estimated total ablation using an empirical relationship between total glacier ablation by Miles et al. (2021) and glacier thinning rates by Hugonnet et al. (2021).
# Calibrating many glaciers takes a lot of time.
# Reduce the amount of time required and only calibrate the glaciers that are
# larger or equal than 0.2 km2. That represents 83% of the glacier volume in the
# Atbashy basin but reduces the number of glaciers to be calibrated by 65%.
# Or, in the cas of an Area threshold of 0.5 km2, the glacier volume covered is
# 59% and the glacier number reduction is 92%
rgi <- rgi |> mutate(Vol_km3 = glacierVolume_RGIF(Area))
rgi_cal <- rgi |>
dplyr::filter(Area >= 0.2)
temperature_cal <- hist_obs |>
mutate(date = as_date(date)) |>
dplyr::select(c(date, starts_with(rgi_cal$RGIId)))
observed_cal <- observed_melt |>
rename(RGIID = RGIId,
totAbl = Qgl_m3a,
year = start) |>
mutate(year = year(year)) |>
dplyr::filter(RGIID %in% rgi_cal$RGIId)
glaciers_cal <- rgi_cal$RGIId
DDFfit <- NULL
for (id in c(1:length(glaciers_cal))) {
cat("Fitting temperature index model for glacier", glaciers_cal[id], "id =",
id, "\n")
DDFtemp <- c(MF = 20, tth = -20)
cat(" Start Parameters: MF =", DDFtemp[1], "Tth =", DDFtemp[2], "\n")
tempGA <- ga(type = "real-valued",
fitness = glacierDischargeRMSE,
temperature = temperature_cal,
observed = observed_cal,
index = id,
lower = c(0.1, -30), upper = c(50, 10),
popSize = 200, maxiter = 100,
parallel = (parallel::detectCores()-6))
cat(" Result: RMSE = ", -tempGA@fitnessValue, "MF =", tempGA@solution[1],
"Tth =", tempGA@solution[2], "\n")
DDFfit <- rbind(DDFfit, tibble(MF = tempGA@solution[1],
Tth = tempGA@solution[2],
rsme = tempGA@fitnessValue) |>
mutate(ID = glaciers_cal[id]))
save(DDFfit,
file = "../../atbashy_glacier_demo_data/GLACIERS/TIparameters.RData")
}Recalculate glacier melt with calibrated parameters.
# Fill in the missing parameters
MFmean <- DDFfit$MF |> mean()
Tthmean <- DDFfit$Tth |> mean()
DDFfit <- DDFfit |>
rename(RGIId = ID) |>
full_join(rgi_elbands |> st_drop_geometry() |> dplyr::select(RGIId, ID),
by = c("RGIId" = "RGIId")) |>
mutate(MF = ifelse(is.na(MF), MFmean, MF),
Tth = ifelse(is.na(Tth), Tthmean, Tth))
MFvec <- DDFfit$MF
names(MFvec) <- DDFfit$ID
Tthvec <- DDFfit$Tth
names(Tthvec) <- DDFfit$ID
melt_cal <- glacierMelt_TI(temperature = hist_obs |>
dplyr::select(-date),
MF = MFvec,
threshold_temperature = Tthvec) |>
mutate(date = hist_obs$date) |>
relocate(date, .before = where(is.numeric))
melt_a_eb <- melt_cal |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt))|>
ungroup()
melt_a <- melt_a_eb |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(Hyear, RGIId) |>
summarise(date = as_date(first(date)),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt)) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |>
left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(),
by = "RGIId") |>
# Assume constant glacier area over the comparison period.
mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6,
.lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),
.lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a),
.ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))
ggplot() +
geom_ribbon(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6,7)]),
aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a,
fill = RGIId), colour = NA, alpha = 0.1) +
geom_ribbon(data = melt_obs_a,
aes(start, Qgl_m3a,
ymin = .lb_Qgl_m3a,
ymax = .ub_Qgl_m3a,
colour = RGIId, linetype = "obs", fill = RGIId),
size = 0.2, alpha = 0.1) +
geom_line(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6, 7)]),
aes(date, Qgl_m3a, colour = RGIId, linetype = "sim")) +
geom_line(data = melt_obs_a, aes(start, Qgl_m3a, colour = RGIId,
linetype = "obs")) +
labs(x = "Date", y = "Melt [m weq/a]") +
scale_linetype_manual(name = "Source",
values = c("sim" = 1, "obs" = 2)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ggtitle(paste0("MF: ", MF, "Tth: ", threshold_temperature)) +
theme_bw()
The calibration is done using an R implementation of a genetic
algorithm Scrucca (2017). The optimization
is slow but it is necessary to use a global optimizer
(nlminb of the stats package for example will not converge
for all glaciers).
Glacier mass balance
The following figure shows the components of the glacier mass balance for a few glaciers in the Atbashy basin.
glacier_balance <- glacierBalance(melt_a_eb = melt_a_eb,
rgi_elbands = rgi_elbands)
glacier_balance <- glacier_balance |>
mutate(.lb = ifelse(Variable == "A_km2",
ifelse(Value*(1-error_stats$sA)>0,
Value*(1-error_stats$sA), 0),
ifelse(Variable == "V_km3",
ifelse(Value*(1-error_stats$sV)>0,
Value*(1-error_stats$sV), 0),
ifelse(Variable == "Q_m3a",
ifelse(Value*(1-error_stats$sQglmelt)>0,
Value*(1-error_stats$sQglmelt), 0),
ifelse(Variable == "Qimb_m3a",
Value*(1-error_stats$sImbal), NA)))),
.ub = ifelse(Variable == "A_km2",
Value*(1+error_stats$sA),
ifelse(Variable == "V_km3",
Value*(1+error_stats$sV),
ifelse(Variable == "Q_m3a",
Value*(1+error_stats$sQglmelt),
ifelse(Variable == "Qimb_m3a",
Value*(1+error_stats$sImbal), NA)))))
ggplot(glacier_balance |>
dplyr::filter(RGIId %in% glaciers[6:7],
Hyear > min(Hyear) & Hyear < max(Hyear),
Variable %in% c("A_km2", "V_km3", "Q_m3a", "Qimb_m3a"))) +
geom_ribbon(aes(Hyear, ymin = .lb, ymax = .ub, fill = RGIId),
alpha = 0.2, colour = NA) +
geom_line(aes(Hyear, Value, colour = RGIId)) +
facet_wrap("Variable", scales = "free_y") +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
theme_bw()
Glacier development and discharge.
We now have glacier discharge (Q_m3s) and the unsustainable contribution to glacier discharge, the imbalance ablation (Qimb_m3a) which is negative for glacier loss and positive for growing glaciers. We are only interested in the contribution of imbalance ablation to river discharge, that is, only the negative part of Qimb_m3a is relevant to us.
Future glacier melt
scenarios <- names(fut_sim)
fut_melt_d <- NULL
fut_melt_a_eb <- NULL
fut_glacier_balance <- NULL
for (scenario in scenarios) {
temp_melt_d <- glacierMelt_TI(temperature = fut_sim[[scenario]] |>
dplyr::select(-Date),
MF = MFvec,
threshold_temperature = Tthvec) |>
mutate(Scenario = scenario,
date = fut_sim[[scenario]]$Date) |>
relocate(date, .before = where(is.numeric)) |>
pivot_longer(-c(date, Scenario), names_to = "ID",
values_to = "Melt")
temp_melt <- temp_melt_d |>
mutate(Hyear = hyear(date)) |>
# Sum melt per year
group_by(Hyear, Scenario, ID) |>
summarise(date = first(date),
Melt = sum(Melt)) |>
ungroup()
fut_melt_d <- rbind(fut_melt_d, temp_melt_d)
fut_melt_a_eb <- rbind(fut_melt_a_eb, temp_melt)
fut_glacier_balance <- rbind(
fut_glacier_balance,
glacierBalance(melt_a_eb = temp_melt,
rgi_elbands = rgi_elbands) |>
mutate(Scenario = scenario))
}
fut_melt_amean <- fut_melt_a_eb |>
# Mean melt rate per glacier
# mutate(RGIId = gsub("\\_\\{d}+$", "", ID)) |>
# group_by(Hyear, Scenario, RGIId) |>
# summarise(Melt = mean(Melt)) |>
# Mean melt rate in the catchment
group_by(Hyear, Scenario) |>
summarise(date = as_date(first(date)),
Melt_mma = mean(Melt)) |>
ungroup() |>
separate(Scenario, into = c("GCM", "SSP"), sep = "_") |>
mutate(SSP = factor(SSP, levels = c("ssp585", "ssp370", "ssp245", "ssp126"))) |>
# dplyr::filter(!str_detect(GCM, "UK")) |>
dplyr::filter(!str_detect(GCM, "IPSL"),
Hyear < max(Hyear)) |>
group_by(Hyear, SSP) |>
mutate(lb = min(Melt_mma),
ub = max(Melt_mma),
mean = mean(Melt_mma)) |>
ungroup()
melt_amean <- melt_a |>
dplyr::select(Hyear, date, Melt) |>
group_by(Hyear) |>
dplyr::summarise(date = first(date),
Melt_mma = mean(Melt)) |>
ungroup()
ggplot(fut_melt_amean) +
geom_ribbon(aes(date, ymin = lb, ymax = ub, fill = SSP), alpha = 0.2) +
geom_line(aes(date, mean, colour = SSP)) +
geom_line(data = melt_amean, aes(date, Melt_mma)) +
# geom_smooth(aes(date, mean, colour = SSP), method = "loess", se = FALSE) +
# geom_smooth(data = melt_amean, aes(date, Melt_mma), method = "loess",
# se = FALSE) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ylab("Melt [mm/a]") +
xlab("Year") +
theme_bw()
Glacier balance over climate scenarios
fut_glacier_balance <- fut_glacier_balance |>
dplyr::filter(!str_detect(Scenario, "IPSL")) |>
mutate(GCM = gsub("\\_ssp\\d{3}$", "", Scenario),
SSP = str_extract(Scenario, "ssp\\d{3}$")) |>
group_by(Hyear, RGIId, SSP, Variable) |>
mutate(.lb = min(Value, na.rm = TRUE),
.ub = max(Value, na.rm = TRUE),
.m = mean(Value, na.rm = TRUE),
# deal with missing values
.lb = ifelse(is.infinite(.lb), NA, .lb),
.ub = ifelse(is.infinite(.ub), NA, .ub),
.m = ifelse(is.infinite(.m), NA, .m)) |>
ungroup()
fut_glacier_balance_basin <- fut_glacier_balance |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear),
Variable %in% c("A_km2", "V_km3", "Q_m3a", "Qimb_m3a")) |>
group_by(Hyear, Variable, Scenario) |>
summarise(Value = sum(Value, na.rm = TRUE),
GCM = first(GCM),
SSP = first(SSP)) |>
mutate(SSP = factor(SSP,
levels = c("ssp585", "ssp370", "ssp245", "ssp126"))) |>
group_by(Hyear, Variable, SSP) |>
mutate(.lb = min(Value, na.rm = TRUE),
.ub = max(Value, na.rm = TRUE),
.m = mean(Value, na.rm = TRUE),
# deal with missing values
.lb = ifelse(is.infinite(.lb), NA, .lb),
.ub = ifelse(is.infinite(.ub), NA, .ub),
.m = ifelse(is.infinite(.m), NA, .m)) |>
ungroup()
ggplot() +
geom_ribbon(data = fut_glacier_balance_basin,
aes(Hyear, ymin = .lb, ymax = .ub, fill = SSP),
alpha = 0.2, colour = NA) +
geom_line(data = fut_glacier_balance_basin,
aes(Hyear, .m, colour = SSP)) +
facet_wrap("Variable", scales = "free_y") +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
theme_bw()
Uncertainties
The shaded area in the above figure indicates the spread of the GCM ensemble used. The width of the uncertainty band of the model method is at least twice the width of the above uncertainty bands, such that it is near impossible to differentiate between the impact of the different socio-economic scenarios. Dramatic glacier volume reduction cannot be avoided even with the least emissions scenario ssp126.
Peak water
ggplot(fut_glacier_balance_basin |>
dplyr::filter(Variable == "Q_m3a")) +
geom_ribbon(aes(Hyear, ymin = .lb, ymax = .ub, fill = SSP),
alpha = 0.2, colour = NA) +
#geom_line(aes(Hyear, .m, colour = SSP)) +
geom_smooth(aes(Hyear, .m, colour = SSP), linetype = 2, size = 0.4,
method = "loess", span = 0.6, se = FALSE) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ylab("Q [m3/a]") +
theme_bw()
Estimated development of glacier discharge.
According to the model, peak water is already reached for the lesser emissions scenarios. We appear to be approximately at peak water for the higher emission scenarios (for which we are on track).
Glacier volume loss
ggplot(fut_glacier_balance_basin |>
dplyr::select(-GCM) |>
dplyr::filter(Variable == "V_km3") |>
dplyr::select(-c(Variable, Scenario, Value)) |>
distinct() |>
mutate(firstm = ifelse(Hyear == 2013, .m, NA),
firstlb = ifelse(Hyear == 2013, .lb, NA),
firstub = ifelse(Hyear == 2013, .ub, NA)) |>
group_by(SSP) |>
fill(firstm, .direction = "down") |>
fill(firstlb, .direction = "down") |>
fill(firstub, .direction = "down") |>
group_by(Hyear, SSP) |>
mutate(V = .m/firstm*100,
.lv = .lb/firstlb*100,
.uv = .ub/firstub*100) |>
ungroup() |>
mutate(SSP = factor(SSP,
levels = c("ssp126", "ssp245", "ssp370",
"ssp585")))) +
geom_ribbon(aes(Hyear, ymin = .lv, ymax = .uv, fill = SSP),
alpha = 0.2) +
geom_line(aes(Hyear, V, colour = SSP)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ylab("V [%]") +
ylim(c(0, 100)) +
theme_bw()
Estimated development of glacier volume.
The reduction of the glacier water reservoirs in the Atbashy basin will amount to about 68%-76% of the current glacier volume.
Next steps
Now we have glacier discharge and glacier imbalance ablation. How do we include these in our modelling workflow? Glacier vignette 4 explains how glaciers are treated in the hydraulic/hydrologic model RS Minerve and details a workflow to implement glacier discharge in a model in RS Minerve.
Open Global Glacier Model (OGGM)
The Open Global Glacier Model is a community project which provides open-source glacier modelling tools in python and educational resources. The modeling framework includes comprehensive tutorials and data-pre and post-processing tools. It can be used as an alternative to the glacier modelling work flow presented here.